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ABSTRACT 

We present the first tests of a new method, the Correlated Component Analysis (CCA) based 
on second-order statistics, to estimate the mixing matrix, a key ingredient to separate astro- 
physical foregrounds superimposed to the Cosmic Microwave Background (CMB). In the 
present application, the mixing matrix is parameterized in terms of the spectral indices of 
Galactic synchrotron and thermal dust emissions, while the free-free spectral index is pre- 
scribed by basic physics, and is thus assumed to be known. We consider simulated observa- 
tions of the microwave sky with angular resolution and white stationary noise at the nominal 
levels for the PLANCK satellite, and realistic foreground emissions, with a position dependent 
synchrotron spectral index. We work with two sets of PLANCK frequency channels: the low 
frequency set, from 30 to 143 GHz, complemented with the Haslam 408 MHz map, and the 
high frequency set, from 217 to 545 GHz. The concentration of intense free-free emission on 
the Galactic plane introduces a steep dependence of the spectral index of the global Galactic 
emission with Galactic latitude, close to the Galactic equator This feature makes difficult for 
the CCA to recover the synchrotron spectral index in this region, given the limited angular res- 
olution of Planck, especially at low frequencies. A cut of a narrow strip around the Galactic 
equator (I &| < 3°), however, allows us to overcome this problem. We show that, once this strip 
is removed, the CCA allows an effective foreground subtraction, with residual uncertainties 
inducing a minor contribution to errors on the recovered CMB power spectrum. 

Key words: methods: data analysis - techniques: image processing - cosmic microwave back- 
ground. 



1 INTRODUCTION 

The Cosmic Microwave Background (CMB) is by far the most 
powerful cosmological probe. The power spectra of its tem- 
perature and polarization anisotropies encode detailed informa- 
tion on the key cosmological parameters. Tremendous experi- 
mental efforts, and especially the currently flying Wilkinson Mi- 
crowave Anisotropy Probe (WMAP, Bennett et al. 2003a, Spergel 
et al. 2006) and the forthcoming PLANCK satellite (Tauber 2004; 
Lamarre et al. 2003; Mandolesi, Morgante, & Villa 2003) will sub- 
stantially advance the sensitivity and resolution of maps of the mi- 
crowave sky. Correspondingly, efficient methods to extract all the 
available information from the measured signal need to be imple- 
mented, exploiting the advances of the data analysis science. One 
of the most challenging tasks is the separation of the astrophysi- 
cal components superimposed on the CMB, usually denominated 
"foregrounds". 

On angular scales larger than about 30' the dominant fore- 
grounds in the relevant spectral region are diffuse emissions from 



our own Galaxy (De Zotti et al. 1999). Synchrotron (Haslam et al. 
1982) and free-free (Haffner, Reynolds, & Tufte 1999, Finkbeiner 
2003) emissions dominate below ~ 60 — 80 GHz (Bennett et 
al. 2003b), while at higher frequencies thermal dust (Schlegel, 
Finkbeiner, Davies 1998; Finkbeiner, Schlegel, Davies 1999) takes 
over. On smaller angular scales foreground fluctuations are dom- 
inated by several populations of extra-galactic sources, with dif- 
ferent spectral behaviour: radio sources, dusty galaxies and the 
Sunyaev-Zel'dovich effect from galaxy clusters. 

A big deal of work has been recently dedicated to develop 
algorithms performing component separation, based on different 
ideas and techniques from signal processing science. Several algo- 
rithms, referred to as "non-blind", assume a perfect knowledge of 
the frequency dependence of sources. The most widely used tech- 
niques exploiting this approach are Wiener Filtering (WF; Tegmark 
& Efstathiou 1996; Bouchet, Prunet, & Sethi 1999) and the Max- 
imum Entropy Method (MEM; Hobson et al. 1998; Barreiro et 
al. 2004; Maisinger et al. 2004; Stolyarov et al. 2005). On the 
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other hand, the emission spectra are generally poorly known. This 
fact has motivated "blind" approaches, which do not make any as- 
sumption on the spectral shape. Particularly promising are meth- 
ods exploiting Independent Component Analysis (ICA) techniques 
(Amari & Chichocki 1998; Hyvarinen 1999), relying on the sta- 
tistical independence of the different components. These methods 
proved to be very effective in extracting the CMB, which is inde- 
pendent of all other sources (Baccigalupi et al. 2000, 2004; Maino 
et al. 2002, 2003; Stivoli et al. 2006; Patanchon et al. 2005). How- 
ever, Galactic emissions are tightly correlated to each other, so that 
the assumption of mutual independence breaks down. Still follow- 
ing the ICA approach, Belouchrani et al. (1997) dropped statistical 
independence for mere uncorrelation, and recovered the missing 
information by exploiting the physically plausible assumption of a 
significant spatial autocorrelation of the individual sources. Bedini 
et al. (2005) also exploited second-order statistics and spatial auto- 
correlation but, at variance with Belouchrani et al. (1997), they gave 
up the assumption of mutually uncorrelated sources. By exploiting 
a parametric knowledge of the mixing matrix, they succeeded in es- 
timating both the mixing matrix and the relevant cross-covariances. 

Here we present the first tests on this method, referred to as 
correlated component analysis (CCA), on a data set as realistic as 
possible. Also, we work either on the whole sky or on spherical 
sky patches, rather than on small plane patches as in Bedini et al. 
(2005). 

The paper is organized as follows. In Section 2 and 3 we de- 
scribe the basic aspects of CCA, and its implementation. In Section 
4 we illustate the tests on simulated skies. In Section 5 and 6 we 
estimate the precision of the mixing matrix and CMB power spec- 
trum recovery. In Section 7 we draw our conclusions. 



2 THE CORRELATED COMPONENT ANALYSIS (CCA) 

We assume that the observed sky radiation is the superposition of 
A'^ different physical processes whose spatial pattern is independent 
of frequency spectrum: 



(1) 



If we have a set of equal resolution observations at M different 
frequencies, the observed signal can be modelled as: 



X = Hs + n, 



(2) 



where x is the A/-vector of observations, H is a M x A'^ mixing 
matrix, s is the A'^-vector of sources and n the M-vector of instru- 
mental noise. The generic element of the mixing matrix is related 
to the source spectra, fc{i^), and to the instrumental frequency re- 
sponse function, bdiy)'- 



hdc = / fc{v)hd{i')dv . 



(3) 



Assuming that the source spectra are constant within the passbands, 
eq. becomes: 



hdc = fc{i^d) / bd{i')du. 



(4) 



meaning that the element hdc is proportional to the spectrum of the 
c-th source at the central frequency Vd of the d-th channel. In any 
case, hdc will be proportional to fc at an effective frequency v^s 
within the d-th sensor passband, depending on both the spectrum 
and the frequency response (see Eriksen et al. 2006). 



If both the mixing matrix H and the source vector s are un- 
known, the problem is unsolvable without additional hypotheses. 
The CCA method exploits information on the second order statis- 
tics of the data. 

The covariance matrix of a generic signal X, defined in a two 
dimensional space with coordinates rf), is: 

C.(r, V) = ([X(e, n) - m] [X(? + r, + V) - m]^), 0) 

where (...) denotes expectation under the appropriate joint proba- 
bility distribution, ji is the mean vector and the superscript T means 
transposition. Every covariance matrix is characterized by the shift 
pair (r, x})), where r and x}) are increments in the f and rj coordi- 
nates. 

From eq. ^2} we can easily derive a relation between the data 
covariance matrix Cx at a certain lag, the source covariance matrix 
Ca at the same lag, the mixing matrix H, and the noise covariance 
matrix Cr,: 



(6) 



This relation allows us to estimate the mixing operator H from the 
covariance matrix Cx of the data. Note that it implicitly assumes 
that the source and the noise processes are stationary. This is true 
all across the sky for the CMB, but only within small sky patches 
for the foregrounds. For this reason, it is convenient to section the 
sky into patches within which foregrounds have approximately uni- 
form properties, and apply the method to individual patches. In this 
way the only remaining issue is the non-stationarity of noise, de- 
pending on the particular scanning strategy adopted, which we do 
not consider in this work. 

If the noise process can be assumed to be signal-independent, 
white and zero-mean, for (r, ijj) = (0, 0) Cn is a diagonal matrix 
whose elements are the noise variances in the frequency channels 
of the instrument, while for (r, ip) / (0, 0) Cn is the null M x M 
matrix. Anyway, if Cn deviates significantly from this ideal model, 
various methods are available to estimate the noise covariance func- 
tion: for example it can be empirically determined using noise maps 
from Monte Carlo simulations. 

Once we have a model for Cn, we only have to calculate Cx 
for a large enough number of nonzero shift pairs (r, x})) to estimate 
both H and Cs. In practice, however, we need to parameterize the 
mixing matrix to reduce the number of unknowns. We note that for 
the scaling ambiguity, we can normalize our mixing matrix to have 
all elements of a reference row equal to 1 . 

The main product of CCA is then an estimate of the mixing 
matrix. Hence this can be considered a "model learning" algorithm. 
Once the mixing matrix have been recovered, source separation can 
be performed with traditional non-blind methods, such as WF or 
MEM, or other Bayesian inversion techniques. 

2.1 Parameterization of the mixing matrix 

To choose a suitable parameterization for H we use the fact that its 
elements are proportional to the spectra of astrophysical sources, 
of which we have some knowledge coming from the theory or 
from complementary observations. The main diffuse components 
present in the PLANCK channels are, in addition to the CMB, the 
Galactic dust, synchrotron and free-free emissions. The frequency 
dependencies of the CMB and of the free-free are known; in terms 
of the anteima temperature we have: 



Ta,cmb(j') oc 



(exp(/i!//fcTcMB) - 1)2 



(7) 
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(8) 



where h is the Planck constant, k the Boltzmann constant, Tcmb = 
2.726 K. Conversely, the frequency scalings of synchrotron and 
dust are not known a priori. The simplest model compatible with 
observations involves only one parameter for each source, the spec- 
tral index fis or the emissivity index (3d, respectively: 

TA,aynch(l') OC , 



(9) 



rA,dust(!^) tX 



exp(/ii//fcrdust) - 1 ' 



(10) 



The mixing matrix accounting for all these sources is then of four 
columns and of as many rows as the number of channels. We are 
able to parameterize the matrix H by exploiting eqs. J3} or J4}. The 
integrals in eq. ^3} can be evaluated with fc replaced by one of the 
emission spectra in eqs. fTl- llOL 

Under the above assumptions, we only have two free parame- 
ters, /3s and /3d- When we work with the high frequency PLANCK 
channels, synchrotron and free-free can be neglected, so the mixing 
matrix has two columns, one for the CMB and one for the dust, and 
only one parameter, l3d- 

Since f3s and I3d vary across the sky, we will apply the CCA 
to sky patches small enough for these indices can be assumed to be 
approximately constant, as we specify in Section 5. 



3 IMPLEMENTATION 

3.1 Computation of the data covariance matrix 

The data covariance matrix defined by eq. ^5) is the fundamental 
tool used to identify the mixing operator. If our data are sampled 
in P pixels, labelled with coordinates rj), we can compute an 
estimate of Cx as: 



Cx(r,V) = 



(11) 



The shift pair (r, ^p) defines a vector that links each pixel to a 
shifted one. In the case of a uniformly sampled rectangular sky 
patch we have a very easy way to define the shift pairs (r, i/;). We 
choose 77) as cartesian coordinates whose axes are parallel to the 
sides of the rectangle, and define a collection of Np shifts, each of 
p pixels, labelled by an index n(n — 0, 1, ..Np — 1): {t„ = npu^} 
and {tpn = npu,,}. From all the combinations of Tn and we get 
Np X Np shift pairs. 

The data we are working with are very different from this sim- 
ple case: they are patches extracted from Healpix (Gorski et al. 
2005) all sky maps, so they are sampled in a sphere rather than 
in a plane and also the grid is not regular. We note that even in the 
simplest case the selected region is not exactly rectangular because 
of the diamond shape of pixels and of the surface curvature. 

To apply the method we need to map the selected patch into 
a geometrically identical one, shifted by (r, tp). This can be done 
only in an approximate way with the Healpix pixelization, and is 
easier in the equatorial region. In the present application we study 
Galactic foregrounds where they are more intense, i.e. at not too 
high Galactic latitudes. We therefore use Galactic coordinates and 
refer to the grid defined by Galactic parallels and meridians to cal- 
culate shifts. 

With the Healpix ring ordering scheme, pixels with subse- 
quent indices are subsequent in longitude, so that, for any inte- 
ger p, a set of pixels imin ^ i ^ imax simply map into the set 



imin + p ^ i ^ imax + p, shifted by p pixels in the longitudinal 
direction. The shift in latitude is trickier. We proceeded associating 
to each pixel the shifted one closest to having the same longitude 
and the latitude increased by Ab = p ■ ds, where ds is the mean 
pixel size. Clearly, in this case A6 is not exactly equal for all shifted 
pixels. However, our shifts are rather small (see below), so that, if 
we are not too close to the Galactic poles, the approximation is 
sufficiently good. 

To choose convenient values for the number of shifts Np and 
the step p, we should care both about efficiency and conditioning. It 
is obvious that, if A'p is large, the method becomes computationally 
demanding. On the other hand, a too small A'p can make the prob- 
lem ill conditioned, thus leading to a lack of convergence. The step 
p should be chosen in order to avoid both redundancy (if p is small, 
some of the covariance matrices are nearly equal) and degeneracy 
(if p is too large, some covariance matrices vanish). In practice, 
the choice can be made empirically, for any p, by increasing A'p 
progressively until convergence is reached. 



3.2 Tlie minimization procedure 

To solve eq. ^6} and estimate the parameters that identify the mixing 
matrix and the source covariance matrices, we minimize the resid- 
ual between the theoretical quantities based on the proposed solu- 
tion and the corresponding quantities evaluated empirically from 
the available data. Our solution is given by: 

(r,E(:,:))=argmin^ || H(r)C.[E(r, V)]H^(r) - 



-Cx(r,v^) + C„(r,,/>) II, 



(12) 



where F is the vector of all parameters defining H, and E(:, :) is 
the vector containing all the unknown elements of the matrices Cs 
for every shift pair. 

To perform the minimization, we used simulated annealing 
(SA, Aarts & Korst, 1989), which exploits an analogy between the 
way in which a metal cools and freezes into a minimum energy 
crystalline structure and the search for a minimum in a more gen- 
eral system. 

The major advantage of SA over other methods is its ability to 
avoid becoming trapped at local minima, which can be very nasty 
in our case. The algorithm employs a random search which not 
only accepts changes that decrease the objective function /, but 
also some changes that increase it. The latter are accepted with a 
probability: 



p = exp(-5//r), 



(13) 



where Sf is an increment in / and T is a control parameter, known 
as the system "temperature". 

Avoidance of local minima is anyway dependent on the "an- 
nealing schedule": the choice of the initial temperature, how many 
iterations are performed at each temperature, and how much the 
temperature is decreased at each step as cooling proceeds. 



4 TESTS ON SIMULATED SKIES 

For this set of tests we used the specifications of the PLANCK mis- 
sion (see Table0. The simulated sky contains: a) synchrotron emis- 
sion as modelled by Giardino et al. (2002), allowing for a spatially 
varying spectral index; b) thermal dust emission (Finkbeiner et al. 
1999), with dust at two temperatures and two emissivity indices; 
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Table 1. Planck specifications. We assume spatially uniform Gaussian noise at the mean level expected for the nominal mission (14 months) 



Center frequency (GHz) 


30 


44 


70 


100 143 


217 


353 


545 


857 


Angular resolution (arcmin) 


33 


24 


14 


9.5 7.1 


5.0 


5.0 


5.0 


5.0 


Rms pixel noise AT (/xK thermodynamic) 


5.5 


7.4 


12.8 


6.8 6.0 


13.1 


40.1 


401 


18291 



c) free-free traced by the Ha emission (Dickinson et al. 2003) and 
corrected for dust absorption with the 100 /im maps from Schlegel 
et al. (1998); d) a CMB Gaussian realization corresponding to the 
best fit WMAP theoretical power spectrum from first year data. We 
produced 100 sets of Monte Carlo simulated maps in the PLANCK 
channels, with different realizations of the CMB and of Gaussian 
noise. 

To test the CCA ability to recover the spectral parameters of 
foregrounds, we used the PLANCK channels from 30 to 143 GHz. 
The resolution of all the maps had to be degraded to that of the 30 
GHz channel (33')- To test the performances achievable with the 
full Planck resolution of 5', we repeated the analysis with the 
high frequency channels from 217 to 545 GHz, where the sources 
to account for are only CMB and dust. 

The noise maps we initially added to each channel were Gaus- 
sian (with the rms levels reported in TablQ and uncorrelated: in 
this case, the only noise term in eq. j6j is Cn(0, 0), which is a diag- 
onal matrix whose elements are the noise variances for each chan- 
nel. The smoothing process applied to degrade the channels to the 
30 GHz resolution not only changes the noise variances, but also 
introduces noise correlation, so in principle the terms Cn(T, t/j) do 
not vanish for (r, ?/)) 7^ (0,0). Nevertheless, we carried out our 
tests assuming uncorrelated Gaussian noise: we then consider only 
the diagonal matrix Cn(0,0), whose elements are the variances 
measured after smoothing each noise map to the 33' resolution. 
This approximation is not the best we could do in dealing with 
noise, but it was purposely adopted to investigate the effect of er- 
rors on noise modelling. 



4.1 The low frequency channels (LF) set 

The LF set includes 5 PLANCK channels, centered at 30, 44, 70, 
100, and 143 GHz. Since we have 4 sources (CMB, synchrotron, 
free-free, and thermal dust) the mixing matrix has 4 columns and 
we want to recover the synchrotron and the dust spectral indices. 
Since all the maps have been degraded to a resolution of 33', we 
operate with the Healpix parameter NSIDE=512, corresponding to 
a pixel size of about 7'. 

The number of shifts allowing a good conditioning of the 
problem is found to be Np — 5. The value of Np and the pixel 
size constrain the amplitude of the shifts that can be used to cal- 
culate the covariance matrices: the minimum shift must correspond 
to one pixel (~ 0.1°); the maximum one cannot generally exceed 
1° to ensure that all covariance matrices are non-null. We chose 
p = 4, which corresponds roughly to our beam size. 

There is a second effect of the resolution: to have sufficient 
statistics, the number of pixels per patch has to be at least ~ 10^. 
With the adopted pixel size, this corresponds to a patch area of 
1500 deg^, which is not optimal for the reconstruction of the spec- 
tral indices that may vary widely with the position. In the syn- 
chrotron template we use here (Giardino et al. 2002), 10% varia- 
tions occur on scales of about 10 degrees. 

We used patches of {Al,Ab) = (50°, 30°) for the lowest 



latitudes and increased the longitudinal dimension for higher lati- 
tudes to roughly preserve the patch area. We sectioned the sky with 
patches centered at longitudes Ic = {0° , 40° , 80° , . . .320° } and in- 
creasing latitudes. In total we analyzed the ~ 80% of the sky with 
|6| < 55°. 

From the first tests we learned that to get a good spectral in- 
dex reconstruction we need a broad frequency range. This can be 
achieved by taking into account additional information from other 
surveys. To this end, we included in our analysis the 857 GHz 
Planck map and the Haslam et al. (1982) 408 MHz map, taken 
as dust and synchrotron templates, setting to all the elements of 
the mixing matrix except the synchrotron one at 408 MHz, and the 
dust one at 857 GHz. 

4.2 The high frequency channels (HF) set 

In the Planck 217, 353 and 545 GHz channels we can ne- 
glect the synchrotron and free-free emissions, so that we are left 
with CMB and dust, and the only parameter to recover is the 
dust emissivity index Pd- These channels allow us to work at 
the best PLANCK resolution: the beam is of 5', so we can use 
NSIDE=2048, corresponding to a pixel size of about r.7. We 
dissected the sky into patches of size (A/,Afo) = (20°, 20°), 
centered at longitudes Ic = {0°, 20° , 40°, ...340°} and latitudes 
fee = {0°, ±10°, ±20°, ±40°, ±60°}. In total, we then analyzed 
- 85% of the sky. 

In this case, we did not need to include other channels to help 
reconstructing the mixing matrix. We used a shift step p = 5, cor- 
responding to an angular size of 0.14°, and Np = 3. 



5 ESTIMATION OF THE MIXING MATRIX 
5.1 Error estimates 

Since, on purpose, the parameters to be determined do not directly 
reflect those defining the sky model, the errors on our estimates can- 
not be simply derived comparing our results with the "true" values, 
simply because the latter in general do not exist. In particular, the 
synchrotron spectral index is varying in the sky, and the dust emis- 
sion is modeled by a two-component spectrum, so they need to be 
treated as we describe below. 

The output of our method is the estimated mixing matrix nor- 
malized at a reference frequency uo- The elements of the mixing 
matrix corresponding to synchrotron and dust at frequency Vd are: 

hout{syn,Ud) = ( — ) , (14) 

/iout(dust,I/d = — 77 JTTf X 7, (15) 

\uoJ exp(hud/kTdust) - 1 

which were derived from eqs. 01. and (ToJ by assuming 
fed(i^) = S{i' — I'd)- These quantities correspond to the mean ratios 
of synchrotron and dust intensities at frequencies and Vd, within 
the considered patch in the simulated sky: 
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hin{syn,Vd) = f — ) 

Vd^^^^+l exp(/lt/o/fc^dust) - 1 
exp(/i!/d/fcrdust) - 1 ' 



/lin(dust, i/d) = f — ) 



(16) 
(17) 



The "observed" indices and /S^ are now directly comparable 
with the derived indices /3s and /3d. We then define our errors as: 



A/3s = 
A/3d = 



log(/iout(syn, i/d)//iin(syn, i/d)) 



log(i^d/i^o) 
log(/iout(dust, i^d)/^in(dust, ua) 



log(i'd/i'o) 



(18) 
(19) 



5.2 Results for the mixing matrix identification 

FigureQrefers to the analysis of the LF set and shows the mean val- 
ues of A/3s and A/3d, over 100 simulations of the sky, as a function 
of longitude for latitudes b — {0°, ±10°, ±15°}. The error bars, 
generally very small and barely visible in the figure, are the stan- 
dard deviations of A/3s and A/3d from the mean. Their small values 
imply that the errors in the spectral index recovery are mainly sys- 
tematic; for both spectral indices, part of the systematic error comes 
from the fact that the models we assume for the analysis differ from 
those defining the sky model. We clearly have problems with the es- 
timation of the synchrotron spectral index, while the errors on /3d 
are generally small. 

We have checked that this effect is due to confusion between 
synchrotron and free-free. The free-free emission is highly concen- 
trated on the Galactic plane but depresses the mean spectral index 
of the combined emission over entire patches, thus biasing the es- 
timate of the synchrotron spectral index. A better recovery of the 
latter would be possible with high spatial resolution maps, so that 
we have a sufficient number of pixels within a much smaller range 
of Galactic latitudes. 

With the angular resolution adopted here, this effect can be 
minimized by cutting out a strip of a few degrees around the Galac- 
tic plane. The situation is substantially better even with a cut of ±1° 
around the Galactic equator, and improves further if we enlarge the 
cut to ±3° (Fig.|2j. All the tests described below are performed by 
applying this cut, even though its effects are modest when working 
with the HF channels. 

Figure|3|summarizes the results obtained over all the analyzed 
region of the sky with both the LF and HF sets. We computed the 
mean errors on synchrotron and dust indices over 100 Monte Carlo 
iterations for each patch, and the mean and standard deviation of 
errors obtained for patches at the same latitude. The bars are the 
standard deviations around the mean over longitudes. 

For the LF set, while the dust index is always reconstructed 
with an error of ~ 0.05, the error on the synchrotron spectral index 
increases with latitude. This is due to the fact that the substantial 
spatial variability of the synchrotron spectral index, implied by the 
Giardino et al. (2002) model, is increasingly difficult to recover as 
the synchrotron signal weakens with increasing Galactic latitude. 

Thanks to the higher angular resolution of the HF set and to 
the lack of other relevant diffuse components besides CMB and 
dust in this frequency range, the CCA allows us to reconstruct the 
dust spectral index with an error ~ 0.03 over the full latitude range 
analyzed. 



6 ESTIMATION OF ERRORS IN THE CMB POWER 
SPECTRUM 

The CCA procedure can be viewed as the first step (learning) in 
component separation. In a further step, the mixing matrix esti- 
mated by this approach can be input to non-blind separation tec- 
niques. The performances of non-blind techniques are normally 
evaluated as functions of the system noise, assuming perfectly 
known mixing matrices. In this section, we estimate the uncertain- 
ties on the CMB power spectrum induced by the errors in the mix- 
ing matrix resulting from our simulations. Since at the moment the 
optimal component separation method exploiting the information 
recovered by CCA has still to be developed, we try to estimate 
the errors in the CMB spectrum under conservative hypotheses. 
Clearly, the errors will depend on the approach adopted to separate 
the individual components. We investigated the general case where 
component separation is performed by a generic linear filter and, 
in particular, by a Wiener filter, and a pseudo-inverse reconstruc- 
tion. The data model we assume provides a space-invariant mixing 
matrix H, obtained with spectral indices /3^ and /3^ as defined in 
Section l5±1 To evaluate the errors caused by an approximated es- 
timate of the mixing matrix, we assumed the estimated matrix, H, 
to be generated by a distribution of the spectral indices such as the 
one described in the above sections. Following a theoretical deriva- 
tion, the errors are evaluated in terms of quantities that are known 
for our simulation, thus allowing the estimation to be made without 
actually separating the components. 

Given the mixing matrix, one way to solve eq. J2j is to work 
in a conjugate space. If we work with sky patches that are small 
enough for curvature effects can be neglected, the correct basis 
functions are the Fourier complex exponentials; on the whole ce- 
lestial sphere, the good basis functions are the spherical harmonics; 
finally, in the presence of an incomplete sky coverage, pixelization, 
and position-dependent noise, the basis functions can be calculated 
according to Tegmark (1996). For simplicity, we are going to derive 
the equations in spherical harmonics. 

In the harmonic space, the problem stated in eq. J2j simply 
becomes: 



(20) 



where the vectors xj^, sjm, and n;^ contain the harmonic coef- 
ficients of channels, sources and instrumental noise, respectively. 
Using a linear approach to component separation, an estimate of 
the vector Si,„, say sim, can be obtained as: 



(21) 



Sim — ^Im ^ 

where W*^-' referred to as the reconstruction matrix, is some 

H 

matrix-valued function depending on the estimate H of the mix- 
ing matrix. The random matrix H is estimated with the method 
proposed in this paper, and is a function of the estimated spectral 
indices. Following Tegmark & Efstathiou (1996) and Bouchet & 
Gispert (1999), the estimation error on the source spectra, in ma- 
trix notation, can be expressed as: 



AC*-') = ^— V 
H 21 + 1 ^ 

m — — l 



(22) 



where (...) indicates expectation and f indicates conjugate transpo- 
sition. By using eqs. <20> and <2H . after a straightforward calcula- 
tion, we have: 

AC<^) = (w|^>H-I)Cs*')(W^)H-I)t± 
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Figure 1. Absolute errors in the spectral index reconstruction for the LF channels without cuts 



+ W^>C„('Hw|V, (23) 

where I is the identity matrix, and Cs and Cn are, respectively, the 
source and the noise power spectra. 

can be simplified by assuming that the source and the 
noise processes are mutually uncorrelated: in this case Cs, C„ and 
AC^^ are diagonal matrices. If we define: 

C^Pii) = C,(')(i,i), (24) 
= Cj\i,i), (25) 

the estimation error on the i-th source is: 

JV 

AC7f(i) ^ AC;f'(z,i) = ^|(w|f)H-I),,fc("(j) + 

M 

where A*' is the number of sources and M the nuinber of channels. 
The error ^C^-\i) refers to the frequency at which the mixing 
matrix has all its elements equal to 1. Eq. y6j highlights how the 
different terms affect the estimation error: the first term accounts 
for the contamination due to the other sources, while the second 
one accounts for the effect of the instrumental noise. Note that in 
the first term of eq. <26> there is the product of the reconstruction 
matrix depending on the estimated mixing matrix H, with 

the true mixing matrix H. As we will see, this term is minimum 
when H = H and increases as the estimated mixing matrix differs 
form the true one. From eq. <26> . we can approximately evaluate the 



error on the spectrum of any component, provided that the source 
and noise spectra are known, without explicitly reconstructing that 
component. 

Tegmark & Efstathiou (1996) and Bouchet & Gispert (1999) 
performed an accurate analysis of eqs. <23t and <26> when the re- 
construction matrix is a Wiener filter: 

W^^ = Cs*''H^[HCs*'^H^ + Cn*'^]~\ (27) 

where the mixing matrix is assumed as perfectly known, that is, 
H = H. They found that, although the Wiener filter reduces the 
error due to instrumental noise, the contamination error due to the 
other sources is always present, even if the actual matrix H is 
known. Moreover, the spectra estimated from the sources, as given 
by eq. J21> . are biased. The following equation can be easily de- 
rived: 

C^''-^ E(§-§L)=W^'HC.(". (28) 

m— — Z 

Bouchet & Gispert (1999) proposed to use a quality factor derived 
from the quantity W^'h in order to correct the estimated CMB 
power spectrum. 

If we reconstruct the sources by the following pseudo-inverse 
filter: 

Wh = (H^H)~'h^, (29) 

we have WhH — I. Then, eq. <23t tells us that the estimation er- 
ror is only due to instrumental noise. Moreover, the mean of the 
estimated power spectrum is: 
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Figure 2. Absolute errors in the spectral index reconstruction for LF channels cutting out the Galactic plane region \b\ ^ 3° 




Figure 3. Mean absolute spectral index errors versus Galactic latitude for the analysis of LF and HF channels with the Galactic cut of ±3° 



Cs'') = a'') + WhC^') (Wj/)t . (30) 

It is apparent ttiat this estimate is unbiased, but for large multi- 
poles the noise term becomes dominant on the reconstructed source 
spectrum. When the noise spectrum is known, as in the case of the 
Planck experiment, eq. <30t can be used to correct the estimated 
spectra. 

To evaluate the uncertainties on the estimated CMB spectrum 
when the mixing matrix is not known exactly, we used a Monte 
Carlo approach: we generated 100 mixing matrices H drawn from 



the probability distribution for the spectral indices described below. 
For each matrix we computed the AC^' (cmb) from eq. <26> . For 
the LF set, these errors have been evaluated at the reference fre- 
quency of 100 GHz using four sources (CMB, synchrotron, dust 
and free-free) and five channels (30, 44, 70, 100, 143 GHz). For the 
HF set we used two sources (CMB and dust) and three channels 
(217, 353, 545 GHz), with the reference frequency at 217 GHz. For 
each channel set, we used the all-sky power spectra of the reference 
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source templates as the input spectra Cf^ (i). The Ci'' (i) for each 
channel have been calculated as: 

Ci')(i) =47ra?/M, (31) 

where Ni is the number of pixels, and the rms pixel noise values Oi 
have been taken from TableQ 

As we have been working with sky patches, we should per- 
form Monte Carlo simulations patch by patch, using the appropri- 
ate basis, and then combine the information over the whole sky. 
However, since we only want to get an indicative amplitude for 
AC^' (cmb), we only performed our trials on a single set of all-sky 
simulations, decomposed into spherical harmonics. On the other 
hand, since an adequate strategy to combine the results obtained 
on patches to achieve a coordinated all-sky component separation 
has yet to be developed, we will confine ourselves to the multipole 
range constrained by the patch size. 

For each iteration of our Monte Carlo simulations, we gener- 
ated the mixing matrix from synchrotron and dust spectral indices 
that are Gaussian-distributed around their true values and f3^. 
In particular, we used standard deviations of 0.1 and 0.05, respec- 
tively, for the distributions of synchrotron and dust indices in the 
LF channels, and a standard deviation of 0.03 for the dust index in 
the HF channels. These values are upper limits to the ones shown in 
Fig.|3| We then exploited eq. <26t to compute the error on the CMB 
power spectrum for both the Wiener filter and the pseudo-inverse 
filter. From our Monte Carlo simulations, we computed the average 
value of ACj^' (cmb), which is the total error on the CMB power 
spectrum estimation. We also estimated the standard deviation ctac 
of the random variable: 

ACf (cmb) - AC'j!j' (cmb) , (32) 

where AC]^' (cmb) is obtained for H = H, and assumed it as the 
error on the CMB spectrum due to the estimation of the mixing 
matrix. 

The results are shown in Figs. |4| and |5| In Fig. |4| we show 
the errors for both the Wiener filter and the pseudo-inverse filter 
from the LF channel data, as functions of the multipole /, together 
with the original CMB power spectrum. Figure |5|displays the cor- 
responding results for the HF channels. All the errors are computed 
without correction for the noise contribution; cosmic variance is 
not included. Note that the Wiener filters are more accurate, but the 
differences with the pseudo-inverse solutions decrease once correc- 
tions for noise are applied. This is particularly true for the HF set, 
where the noise level is higher. The total error on the CMB power 
spectrum AC^^ (cmb) turns out to be of the order of 1% for both 
the channels sets, and decreases from lower to higher multipoles. 
In both cases the contribution of the errors in the estimation of the 
mixing matrix is subdominant. 



7 CONCLUSIONS 

The tests described here demonstrate that the Correlated Com- 
ponent Analysis (CCA) method (Bedini et al. 2005), applied to 
simulated data with PLANCK specifications, is a promising tool 
to estimate the mixing matrix parameterizing the frequency scal- 
ing of different astrophysical signals in cosmic microwave back- 
ground (CMB) observations. For the low frequency PLANCK chan- 
nels, from 30 to 143 GHz, and with the angular resolution of 
33', the most accurate estimates are obtained in the latitude range 
[—30°, +30°], where the foregrounds we aim to recover are rela- 
tively strong. 
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Figure 4. Estimated errors on the CMB power spectrum at 100 GHz for the 
LF set, for Wiener filter and pseudo-inverse filter. For each filter, the thicker 
(upper) line shows the total error, the lighter (lower) one is the contribution 
due to errors in the mixing matrix estimation. 
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Figure 5. Same as in Fig.|3 but for the HF set, at the reference frequency 
of 217 GHz. 



However, the strong concentration of intense free-free emis- 
sion in a narrow strip (a few degrees wide) around the Galactic 
equator introduces a large gradient of the global spectral index of 
Galactic radio emission (synchrotron plus free-free) that prevents 
an accurate recovery of the synchrotron spectral index with the low 
resolution of the LF channels. The problem is cured cutting out 
the strip at \b\ < 3°. Having done that, both the synchrotron and 
the dust spectral index are recovered with a mean absolute error of 
0.05. At higher latitudes, the mean error on the synchrotron index 
increases to 0.08, while that on the dust index is unchanged. 

The situation is even better with the high frequency PLANCK 
channels, from 217 to 545 GHz, thanks to the better spatial resolu- 
tion (5') and to the fact that the only relevant foreground, at least 
on intermediate angular scales, is Galactic dust. 

We have shown that, exploiting non-optimized component 
separation techniques, such as Wiener filter and pseudo-inverse fil- 
ter, such errors allow us to estimate the CMB power spectrum with 
an uncertainty of the order of 1% on the angular scales constrained 
by our patch size. 
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We conclude that the CCA can be a promising independent 
way to probe the spectral behaviour of the main foregrounds affect- 
ing the CMB. This second order statistical approach may allow us 
to increase our knowledge on foregrounds and to perform compo- 
nent separation with traditional non-blind methods with a minimum 
number of priors. 
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